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1. Introduction 

The dynamics of populations and species in the context of biological or ecological 
systems have attracted considerable attention in the community of statistical physicists 
in the last decade. Though there has been much progress PH E], it is rare that the 
behavior of macroscopic populations can be predicted, even for "simple models," from 
a set of stochastic rules for an individual's propensity to survive and/or reproduce. 
The difficulties can be traced not only to the presence of many degrees of freedom, 
with complex internal interactions (e.g., mutualistic or predator-prey), but also to 
the non-trivial couplings to external reservoirs (such as energy or food). As a result, 
even if we only focus on systems in stationary states, we must recognize that these 
are non- equilibrium steady states, so that the well-known methods of equilibrium 
thermodynamics should not be blindly applied. At present, given the absence of a 
universally applicable framework of non-equilibrium statistical mechanics, progress is 
made through understanding "one system at a time." Within this context, we recently 
studied a model of coevolution H] , based on the one ( "tangled nature" ) introduced 
by Hall, Christensen, and collaborators El Hj • The motivations behind these studies 
are varied. An early model for coevolution and speciation, introduced by Bak and 
Sneppen jS], consists of competing species, according to a preassigned notion of "fitness." 
Speciation arises by imposing a crude version of "mutation," i.e., letting the least fit 
species (as well as some of their "neighbors") be replaced by new species with different, 
randomly chosen fitness. Despite their simplicity, such models appear to settle into 
a steady state which exhibits avalanches of extinctions. Though they are hailed as 
showing a link between Darwin's principle of "survival of the fittest" and Gould's notion 
of "punctuated equilibria" |U E3 HI] , these models are thought to be too simplistic in 
at least two aspects. Firstly, mutation and selection act on individuals of a populations, 
rather than on entire species at once. Secondly, whether a species is "fit" is not a static 
notion, but rather depends on what other species are present in the ecosystem. 

To address both issues, Hall, et.al. 00 Hj recently introduced an individual-based 
model with a dynamically evolving "fitness landscape." The "ecosystem" consists of 
individuals, each of which is said to belong to a "species" identified by a "genome" 
represented by a string of bits. Mutations are built in through the random flipping of 
bits, at a constant slow rate, in the genomes of newborn individuals. In earlier models of 
speciation, an individual's fitness, i.e., its reproduction probability, is purely a function 
of the bit string and fixed for all time. Here, the reproduction probability depends on the 
relative abundance of all other species, so that "fitness" becomes a dynamic concept. 
The "interspecies interactions" are, by contrast, comparatively static in nature and, 
for simplicity, are introduced via fixed quantities associated with pairs of genotypes. 
Thus, the model accounts only for whether (an individual of) a species is beneficial or 
detrimental for another species, regardless of the presence of any other species. Again for 
simplicity, all species reproduce asexually, with identical fecundities. Simulations reveal 
several interesting behaviours, including the presence of long-lived states separated by 
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bursts of high activity - "punctuated equilibria." Unfortunately, this model proved too 
complex for analytic understanding, and only phenomenological descriptions have been 
advanced to date. 

In an effort to gain some insight into how its remarkable properties arise, we 
considered a variation of this model, with simpler interspecies interactions, which 
enables us both to carry out much longer simulations and to perform linear stability 
analysis [3J Hj. In addition to observing a self-similar picture of intermittency or 
"punctuated equilibria" over several decades of generations, we found that standard 
measures of diversity display 1/ / noise in their power spectral densities. This property is 
quite consistent with another finding: The life-time distribution of the long-lived states 
approximately follows a inverse-square power law ^2]- More detailed investigations 
of the very long-lived states, which typically consist of a community of a handful of 
dominant species along with a "cloud of mutants," reveal the reason behind their 
longevity. Essentially none of the closest mutants of the main species in such a 
community are "dangerous," in that their interactions with the parent species inhibit 
their exponential growth. In other words, the original community is (linearly) stable 
against invasion by its most closely related mutants. Only a small fraction of the next- 
closest mutants (with genes differing from a dominant species by two bits) are dangerous, 
leading to the eventual demise of the "quasi-steady" state. Despite the discovery of these 
connections, full analytic understanding is still beyond our grasp. In particular, since 
the theoretical analyses were based entirely on a heuristic, deterministic ("mean-field") 
equation of motion, it is unlikely that they can account for the most intriguing behaviour 
stemming from a stochastic dynamics. 

In the present paper, we address some of the issues associated with a fully stochastic 
description. Starting from a master equation which governs the evolution of all the 
details ("microscopies") of the model, we first demonstrate that the deterministic 
equation in [3] emerges, provided all correlations are ignored. However, it is no easy 
task to find a quantitative understanding describing the quasi-steady states (QSS), in 
which the system appears to be stationary but actually has long, finite lifetimes. One 
complication lies with the inherent metastability aspect of a QSS. Another is that, due to 
mutations, the populations in a QSS consist of two components: a handful of dominant 
species (with at least several hundred individuals, in the specific simulations we ran) 
and a larger number of minor species (with much less than one hundred individuals). 
Our approach to the solution is a two-step process. First, for each QSS, we develop a 
full understanding of a corresponding "truly stationary state" (TSS), i.e., one with only 
the dominant species. Each TSS is associated with an A/"-species fixed point discussed 
in ^] and can be accessed by setting the mutation rate, /i, to zero. The second step is to 
account for O (//) effects, by including a limited "cloud of mutants." Needless to say, a 
careful definition of such a community will be necessary before any analytic progress is 
possible. Thus, we will only take the first step here, deferring the more complex problem 
to a later publication. Beyond that, our eventual goal is to predict the more fascinating 
phenomena, such as power-law distribution of QSS lifetimes, "punctuated equilibria" 
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(or intermittency) , and 1/f noise. 

In the next section, for completeness and the readers' convenience, we briefly review 
the specifications of the model. The following section is devoted to the master equation 
and the derivation of the deterministic evolution equation used in |I] . Fluctuations and 
correlations in a TSS, as well as comparisons to simulations in specific cases, are the 
focus of sections |U and El Conclusions and an outlook can be found in section El The 
Appendix is devoted to some of the technical details. 

2. Model specification and algorithm 

The model we considered jl] is a simplified version of the one introduced in [SI El Hi- 
lt consists of a population of individuals, each of which is associated with a string L 
of bits (0 or 1), representing a "genome" of L "genes." The bit-strings, or genotypes, 
are labeled by integers I , which lie between 1 and A/" max = 2 L . For simplicity, we will 
use the term "species" to distinguish individuals with different bit-strings or genotypes. 
The population evolves asexually in discrete time steps (t = 0,1,...), which may be 
thought of as "years" or "generations." In our model, all individuals of a generation 
die when those of the next generation are "born" (as in, e.g., aphids). Thus, for any 
particular run (or "history"), the system is fully specified by the set of integers ni(t) 
(1 = 1, A/" max ) representing the number of individuals of genotype I in generation t. 
In cases where the explicit index / is not necessary, we will use the "vector" notation: 



To model competition and interspecies interactions, we let an individual die with 
some non- vanishing probability before it reproduces, e.g., salmon that die in the oceans. 
All survivors then give birth to F offspring, which constitute the next generation. 
Competition for resources (e.g., space, energy, food) is often introduced via a Verhulst 
[T5] factor, which also plays the role of preventing unlimited growth and enters typically 
via a ratio N tot (t)/N . Here, A is a parameter representing the "carrying capacity" of 
the "ecosystem," and 



is the total population at time t. With a "healthy" fecundity (F), the population is 
unlikely to "collapse" (A to t = 0). Instead, N tot (t) is rarely far from A . Interspecies 
interactions are modeled by a matrix M, the element Mf being the effect of the species 
J on species /. For reasons provided in jl], we set Mj = and choose random off- 
diagonal elements from a uniform distribution over [—1, 1]. In all our simulations, M 
is fixed at t = and does not evolve in time. If both Mf and Mj are positive, the 
two species are said to be "mutualistic," and if they are of opposite signs, we have 
a predator-prey relationship. Not surprisingly, populations with both elements being 
negative are extremely unstable. Finally, the probability of an individual of species / 



n(t) = {n%(t),...,n Mm ^(t)} . 



(2.1) 




(2.2) 
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to survive to reproduction is specified by jlj 



P(I;n(t)) 



1 



(2.3) 



1 + exp [-M^nj(t)/N tot (t) + N tot (t)/N ] ' 



which is often shortened to just P(I). Here, as in the rest of the paper, we use the 
Einstein summation convention, i.e., a repeated index (e.g., J in M/rij) is summed over 
(from 1 to A/" max , in general). As will be clear below, where we deal with systems with 
large Nq ("macroscopic," though not necessarily in the sense of typical thermodynamic 
systems), quantities with a single subscripted index (e.g., rij) are generally of order 
iVo; those with a single superscript, of order l/N Q ; those with both (e.g., M/) or 
none (e.g., P(I)), of 0(1); etc. Exceptions are noted with a caret (or "hat"). For 
example, we denote the "normalized" covariance matrix for the populations (i.e., 
((riinj) — (rij) (rij)) /No) by Gu, which is a quantity of 0(1) rather than O(Nq). 
Similarly, its inverse (f /J ) is also of 0(1) , as opposed to of 0(1/Nq). To avoid 
confusion, we will remind the readers of such exceptions at the appropriate points. 

The last ingredient in our model is mutation. In the absence of mutations, the 
diversity of the population never increases with time. Indeed, the only rigorously 
stationary state is the collapsed one, rij = for all /. Nevertheless, non-trivial states 
(with extremely long life times and called TSS's here) are typically reached, for all 
relevant time scales. While the rigorously stationary state is independent of the intitial 
condition n(0), the TSS states are entirely dependent on the initial population, and they 
display the dominant characteristics of the corresponding QSS's. Identified by a fixed 
point (FP) of the deterministic evolution equation, 



a TSS provides the basis for linear stability analysis and for future investigations of the 
associated QSS. 

Returning to mutations, they are important, not only to promote diversity and 
to model "speciation," but also to provide the main ingredient for the interesting 
phenomenon of, say, intermittency ("punctuated equilibria") 0]. As in all bit-string 
models involving asexual reproduction, we allow each offspring to carry the genes of 
its parent, except for a probability of fi/L that each bit be changed. As a result, on 
the average and for small /i, a survivor produces fiF offspring with different genetic 
material. To keep track of the "biodiversity," we define the species richness A/"(t) as the 
number of populated species at t (i.e., only M Fs are present at t). Another common 
measure which characterizes the relative abundance of the species better (but will not 
be the focus here) is the Shannon- Wiener index, Ylj P(J) l n P(^); where 



is just the fraction of species J in the system. 

Let us briefly summarize the algorithm used to simulate this model, referring to 
|E] for the details. There are three layers of nested loops: (1) over generations t, (2) 
over Af(t), and (3) over nj{t). In the innermost (last) loop, each individual produces, 



n I (t + l)=n I (t)FP(F ) n(t)}) 



(2.4) 



p(J) = nj(t)/N tot (t) 



(2.5) 
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with probability P(I), F offspring, each of which is allowed to mutate before n(t + 1) 
is recorded. In most of our studies, we chose L = 13, N = 2000, F — 4, /i — 10~ 3 , and 
-^tot (0) = 100, with random initial n(0). Thus, Af max = 8192, though far fewer species 
are typically present in the system at any given time, i.e., J\f <C W max . 



3. Master equation and mean-field theory 

In our recent investigations, the Monte Carlo studies generate stochastic sequences of 
configurations (i.e., non- negative integers, or just "points," in the A/" ma x-dimensional 
space) but the theoretical analysis was based on deterministic, heuristic equations of 
motion for the averages of the populations ("mean- field" theory). To understand the 
full stochastic process, we need a complete description, involving the probability that 
the system is found with a specific number of individuals at time t, namely, V(n,t). 
Its evolution is governed by the master equation, given in an appendix of 4J. Before 
turning to this equation, let us emphasize the difference between V (n, t) here and the 
n(t) above. In the former, n is a co-ordinate (in A/" max -space) and V an evolving 
function in this space. By contrast, n(t) is just the trajectory of a single point in 
this space. A single Monte Carlo run generates a particular trajectory (or "history"): 
ni (t), and can be represented as a (Kronecker) delta function jumping from point to 
point: PaMC run (n, t) = Yli $ { n ii n i (*))■ The full dynamics of V (ft, t) is simulated by 
averaging over many runs, and thus difficult to access. When we turn out attention to the 
TSS's below, the process simplifies, since they are characterized by static distributions: 
V* (n). Then, it is sufficient to perform a single, long run during which the system 
rarely wanders far (say, O (y/Noj) from the neighborhood of a specific point (typically, 
fixed points of the mean- field evolution equations). 

Returning to the issue at hand, finding an equation for V (ft, t), we recapitulate the 
probability for an individual of species / to survive: 



P(I) = \ l + exp 



iV tot M jnj^- ] 



Again, note the different interpretation we give for this expression versus equation ()2.3|) . 
Here, P (I) — P (I; n) denote A/" max functions defined in Af max -space, independent of t. 
In contrast, for a particular MC run at a particular time t, we need only M (t) functions 
(for a specific point {nj(t)}, in a much smaller, AA-space). Next, we must keep track of 
all of the possible number of survivors. Since all of these individuals reproduce, we will 
call them "parents." Defining the symbol [^] as the rate for mj individuals to survive 
from the original nj, we simply write a binomial distribution: 



rii 
mi 



^ — v [p(i)r[i-p(i)r~ mi . (3.2) 

mj! (ni — mj)l 

Next, each parent gives rise to F offspring, not every one of which is of the same 
species. In the simulations, it is possible to have a mutant whose genome differs 
from the parent by two or more bits. However, this is quite rare, being less than 
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O (/J. 2 N tot ) = O (fi 2 N ) ~ 1CT 3 . Thus, we will keep our analysis simple by restricting 
our analysis here to mutations which flip only a single bit. Then, there can be only L+l 
possible varieties of offspring for each parent. To account for these, let us introduce the 
notation 

hjfl for the number of offspring from parent J with no mutations 

6j jCr for the number of offspring from parent J with the ath bit flipped 

and define the multinomial-like symbol 



Fmj 



(Fmj)\ 



^ n 



(3.3) 



This is the probability that the Fmj offspring are distributed into the specific set 
{bj t o,bj t i, ...,bj : i}. The last ingredient needed is the connection matrix 



A J ' a 

K 



1 if genotype K is J with the ath bit nipped 
otherwise 



so that the number of offspring born into species K due to mutations is 



>J,OL 



(3.4) 



(3.5) 



J,a>0 



With these ingredients, we arrive at the master equation p. 

Fmj 



bjfii bjA, bjj 



n 



rii 
mi 



V(n,t) , (3.6) 



V(n',t + 1)= UHn' K ,bK,o + B K )l[ 

n,m,{b} K J 

where 5 (n', n) is the Kronecker delta. Given a particular initial configuration n , V (n, t) 
can be found, in principle, by recursion with V (n, 0) = 5{n,no). In that sense, we 
note that the more precise notation is V (n, t\fio, 0), explicitly showing that it is the 
probability to find our system in state n at time t, conditioned on a specific initial 
condition. However, this notation seems unnecessarily cumbersome, so that we will just 
use V (n, t) in its place. 

Note that equation ()3.6|) is just a special example of the general evolution of 
conditional probabilities in a Markov process, i.e., 

V (n', t + l) = J2R (n'\n) V (n, t) , (3.7) 

n 

where R(n'\n) is the conditional probability for finding the system in ft' given that it 
was in n, also known as the transition rate. For our case, R is explicitly 

Fmj 



R (n'\n) = II 6 ^ b « + B k) II 

m,{b} K J 



lift 



hi 



, -,bj,L 



n 



nj 
mi 



. (3.8) 



| Many master equations are written for continuous time, in the form of dtV(C',t) = 
"YIiqi L (C, C) V(C ', t). For discrete-time processes, we find it more convenient to express the evolution 
in the form used here. 
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Once V (ft, t) is found, the time dependence of the expectation value of any quantity 
can be obtained via 



(•>t = I>Pfa*) 



(3.9) 



in principle. Specifically, our main interest here will be (nr) t , the average number of 
individuals of species I at time t, as well as quantities quadratic in n (e.g., covariances 
and correlations). For example, to see how (nx) t evolves in time, we multiply equation 
()3.7j) by n' K , and sum over n' 

(n K ) t+l = (f K (n)) t , (3.10) 



where 



Ik (n) = n K R $ 



n 



Exploiting the explicit form of R above and 



711 



777 



ml (n — m)\ 



mq 



(3.11) 



(3.12) 



it is straightforward to find /. Reminding the readers of the n dependence in P, we 
write explicitly 

L 

f K (ft) = F n K P (K; n) + (fi/L) ^ ^ IS.fnjP (J; n) . (3.13) 

J a=l 

The interpretation of various terms in this expression is clear: F offspring are born to 
njP (J) survivors of species J, with rearrangements into the new generations due to 
mutations. Inserting it into equation (|3.1()j) . we arrive at an exact equation 

L 

(n K ) t+1 =F {l-y){n K P{K)) t +^/L)Y J Y,^ J K{ n J P i J ))t -(3-14) 

J a=l 

Its simplicity is deceptive, since the P's are, from equation ()3.1j) . non-trivial functions of 
n. If a formal expansion in powers of n were inserted for these functions, then averages 
of all powers, (nji~iK...) t , will appear on the right-hand side. Of course, equations for 
these new averages can be written formally, but the result would be more complex than 
the BBGKY hierarchy ^3]. Needless to say, good approximation schemes are crucial 
for further progress. 

One such scheme, also known as the "mean-field" approximation, is to ignore all 
correlations in order to produce a closed equation for (n K ). Thus, we replace all averages 
of the products by products of the averages: 



so that, e.g. 



(njn K ...) -> (nj) (n K ) 



(n K P{K;n)) -> (n K ) P {K; (n)) . 



(3.15) 
(3.16) 



Fluctuations and correlations in biological coevolution 



9 



In terms of the less cumbersome notation n K (t) = (n K ) t , equation (|3.1U|) reduces to 

n K (t + l) = f K ({nj(t)}) . (3.17) 

Apart from the change of notation (P K ({nj(t)}) — > P(K;n(t))), this equation is 
precisely the starting point of our earlier analysis, equation (2) in jlj: 

n K (t + 1) = n K {t)FP K ({nj(t)})[l — fj] + {^/L)F ^ n I(K) {t)P I{K) {{nj(t)}) , (3.18) 

I(K) 

where / (K) runs over all values of / that differ from K by one bit. 

In the low-mutation-rate (/i <C 1) regime, much of the behaviour of a QSS is well 
approximated by a TSS (/i = 0), which can be understood in terms of the fixed points 
of the mutationless version of this mean-field theory. Here, let us briefly summarize the 
predictions of this theory. With /i = 0, the number of populated species in the system, 
Af, never increases. Therefore, we can restrict our attention to a small subspace (Af- 
dimensional, with Af < 5 in typical simulation runs) of the full 2 L -dimensional space. 
In our previous work [I], we use tildes (e.g., M) to emphasize this aspect. To keep the 
notation simple, we will drop the tildes here and keep in mind that, for example, M is 
an Af x Af matrix. Also, with Af being an 0(1) quantity (i.e., small compared to N ), 
the population of every species is generically O (A ). 

Next, fj simplifies to 

fi(n) ^FnjP(P,n) . (3.19) 

Since the fixed point, n*, obeys 

n* = f(n*) , (3.20) 
we have FP (J; n*) = 1 for all /. Defining the inverse of M by 

W = 1VT 1 , (3.21) 

(with elements Wf) and the sum 

<7 = £V/, (3.22) 
i j 

we found that the fractions of each species are given by 

p*(/) = n*/A t ; t = ^H///a, (3.23) 
J 

and the total N* ot , by 

N:JN = \n(F - 1) + l/a . (3.24) 
Finally, the elements of S, the stability matrix, is 

S J j = f/ (ft) = dfj/dnj^ = 8f + Aj , (3.25) 
where 8j is the unit matrix. Here, 

Aj = (l - (M/ - ln(F - 1) - 2/a) p* (I) . (3.26) 

(denoted by Ajj in |1| and referred to as the community matrix in biological literature) 
plays the role of a "restoring" force, driving the population back towards the (stable) 
fixed point. 
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4. Fluctuations and correlations of populations 

It is clear that all stochastic aspects of the system are lost in the mean-field 
approximation. In this section, we study the simplest aspect, namely, an approximate 
description of the long-lived, quasi-steady states (QSS). Now, as we have seen in 
simulations |lj , these states are dominated by a few mutually supportive species, along 
with a few individuals of closely related "benign mutants." Thus, we restrict ourselves 
in this paper to systems without mutation, so that several simplifications apply. First, 
there is the obvious reduction of R from equation (|3.8j) to 



M m K I 



m 

mj 



(4.1) 



Next, since there can be no new species, we can focus on populations with M ~ 0(1) 
species, each of which having O (N ) individuals. Finally, though it is possible for a 
fluctuation to collapse the entire population (which is, rigorously, the unique stationary 
state associated with equations (j3.7l4.lj0 . such events are so extreme that the life times 
of a typical non-trivial state should be O (e^ ). So, for all practical purposes, we may 
regard such states as "truly" stationary (the TSS's). § 

A complete description of a TSS is provided by a stationary distribution V* (n), 
which satisfies 

V* (n') = J2 R v * (™) • ( 4 - 2 ) 

ri 

Given V* (ft), we can find stationary averages of any quantity: 

(.)* = 5>)P*(n) . (4.3) 



Here, our main interest will be mean populations and their correlations: 

(n/)* and (ninj)* - (raj)* (rij)* . (4.4) 

Even in the absence of mutations and with simplified i?'s, there are substantial 
non-linearities (through P (ft) ) in the problem that prevent us from finding solutions 
to equation (|4.2jl in general. On the other hand, observations in simulations (and 
lessons learned from the central limit theorem) show that our distributions are well 
approximated by Gaussians. Indeed, a systematic expansion for V* can be formulated, 
starting from the Gaussian form: 



1 /2 

V G (ft) = (27rN y Af/2 (det f) exp 



2N ( 



'ni - m) t IJ (rij - nj) 



(4.5) 



where nj and f IJ are parameters to be determined. Of course, the first of these is of 
O (Nq). As for the latter, we expect the fluctuations in our problem to be O (y/No), 

§ It is possible to force such states to be rigorously stationary, by slight modifications of the rates. The 
simplest is to let the survival probabilities P (I) be unity when rij = 1. 
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i.e., covariances of O (No). To incorporate this expectation, we have chosen to write Vq 
in a form such that the matrix T has elements of O (1): 

f 7J ~0(l) , (4.6) 

despite the presence of two superscripts. 

Within the context of such a scheme, it is possible to compute these quantities 
from the microscopic rates. Before proceeding to this computation, we remark that this 
approximation for V* will lead to the predictions 

( ni )* = fii and (n/raj)* - (m)* (rij)* = N G U , (4.7) 

where G is the inverse of T, i.e., 

G u f JK = 8? . (4.8) 

Of course, we expect 

Gu~0(l) , (4.9) 

despite its two subscripts. 

Though such approaches are well known [15J, we provide a few details in the 
Appendix, both for the sake of completeness and for the convenience of readers 
unfamiliar with these methods. Here, we present only the highlights of the analysis and 
how they apply to our case. As shown below, the agreement between our predictions and 
simulation data in three specific cases are excellent and validates this entire approach. 

Turning to the computation of nj and T (or G), we may expect the former to 
be simply related to the n\ of mean-field theory (equation (|3.20j) ). As shown in the 
appendix, we have 

nj = n}[l+O(l/N )} (4.10) 

(provided none of the eigenvalues of S are close to unity). Indeed, the first correction 
can be computed explicitly. Since these corrections are not necessary for finding the 
fluctuations and correlations, we will not quote the results here, but refer the reader to 
equation (jA.22|) in the Appendix. To find G, we need not only // (n) and the stability 
matrix S (given by equations (|3.19J3.25J3.2fi|) ). but also 

H I j(n) = ^2n' I n' J R(n'\n) . (4.11) 



m 

mi 



(4.12) 



In the limit fj, — > 0, R reduces to equation (|4.1|) . so that 
Hu(n) - ( Fm i) ( Fm j) II 

in I 

= F 2 [njnjP (I) P (J) + SjmP (I) (l-P (I))] , (4.13) 

which is indeed of the form in equation (|A.26|) : Hjj = fjfj + N Hjj. Thus, we readily 
identify Hu of equation ()A.30j) : 

H u = 5/F 2 n*P (J; rf ) [l-P (/; rf)] /N Q . (4.14) 
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Figure 1. Theoretical Gaussian (thick, grey curves in the background) and simulated 
(thin, black curves in the foreground) TSS probability densities for Af — 2, F = 4, 
and Nq = 2000, corresponding to the second column in table ^ The Monte Carlo 
simulation was performed over 524,290 generations with zero probability of mutations, 
(a) Contour plot of the joint distribution for n\ and 112, equation 1)4. 5[l . The negative 
slope of the long axis indicates that n\ and n 2 are negatively correlated. (6) Marginal 
distributions for m (right) and n 2 (left), equation l|4.22l) . 



But, FP (J; ft*) = 1 for all I, so that 

H IJ = 5j(F-l)n* I /N Q . (4.15) 

Referring to the Appendix again for the details, G can be obtained from S and H via 
the linear relationship 

G - SGS T = H . (4.16) 

Inserting the explicit form for H into equations (|A.39l IX.40|) . we arrive at a complete 
solution for the covariance matrix in terms of A a , uf, and Vj (respectively, the eigenvalues 
and the left and right eigenvectors of S^, normalized by u^v 1 } = 5^) : 

Gu = (F-1)J2 #XriT + • (4.17) 

a,b,K 

Before comparing these predictions to simulations, we provide an intuitive picture 
for these fluctuations and correlations. Since we are concerned with distributions well 
approximated by Gaussians, the underlying process is just an Ornstein-Uhlenbeck one 
[TB] . Such a process can be described by a Langevin equation: 

e(t + 1) - e(t) = Ae(t) + ff(t) , (4.18) 

where A plays the role of restoring forces which drive e towards zero, and 77 is a Gaussian 
noise with 

(fj(t)} = 0; (ff(t)ff(t')) = H8(t,t') . (4.19) 
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The notation for the two matrices (A, H) is chosen deliberately. For this problem, they 
are precisely those given above: equations ()3.26I4.14|) . while e(t) is just the deviation 
from the average: 

e(t) =n(t) - n . (4.20) 

The deterministic part of equation (|4.18jl is intuitively clear, being the same as in the 
mean-field evolution equation, linearized about the fixed point. Now, fj can be seen 
as the (microscopic) noise in the stochastic process. In our case, this is clearly due to 
the uncertainties associated with survival. Since we have a simple two-state (dying or 
surviving) random process, we can hardly be surprised by the presence of the factors 
nP (1 — P). Also, since this randomness is imposed on each individual, the diagonal 
form of H can be anticipated. Perhaps only the factor F 2 cannot be easily surmised. 

We have carried out simulations for ten cases, with M = 2, 3, and 4. For simplicity, 
in each case we chose a set of species which served as the dominant ones in one of the 
ten QSS communities included in Table I of jlj. Of course, since \x was set to zero in the 
simulations reported here, the communities are TSS's. The runs were carried out for 
524,290 generations each. From the recorded nj (t), we computed the averages (ni) MC 
and the covariance matrix Gmo For the M = 2 case, we can easily display a histogram 
of all the populations (figure da)). In addition, in figure [TJb), we show how good the 
Gaussian approximation is by plotting the projections of both this histogram (onto one 
or the other axis) and the theoretical predictions, e.g., 

^ P roj(ni) = / dn 2 P G (ni,n 2 ) (4.21) 
detf 



27rA r 22 




detr . _ x2 

— (m - nx) 

2A r 22 



(4.22) 



We emphasize that the theoretical curves are produced with no fitting parameters - 
all quantities were computed from the model specifications (i.e., N , F, and M). For 
the Af > 2 cases, it is difficult to display full histograms. Instead, we only provide 
the comparison for the averages and covariance matrices for three particular TSS's in 
t ablest and |U As we see, the agreement is excellent, well within the expected accuracy 
of the approximation, O (I /No), and the statistical errors, O (1/^524,290). 



5. Distributions containing dynamical information 

In the previous section, we focused on the static distribution of the populations in a 
TSS. In other words, we can compute (within the Gaussian approximation) correlation 
functions that involve any number of species, all "at the same time." Here, we investigate 
the information contained in the dynamics of the stochastic process, i.e., time-dependent 
correlations. For an evolving population, a natural question is how the system changes 
from one generation to the next. To probe this issue at a quantitative level, let us 
consider two examples, the statistics of "steps," 

s(t) = H(t+ 1) -n(t) , (5.1) 
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Table 1. Theoretical results for two TSS communities with TV = 2 and 3, respectively, 
compared with corresponding quantities from Monte Carlo simulations with Nq — 2000 
over 524,290 generations. The communities are defined by the interaction matrices 
M. The quantities shown are the normalized fixed-point population vector, n* /Nq 
[equations l|3.23|) and (|3.24|) ] and the corresponding Monte Carlo average (n)Mc/No, 
the normalized population covariance matrix G [defined by equation l|4.7|l and 
calculated from equation l|A.36|) with 20 iterations for TV = 2 and 30 iterations for 
TV = 3] and Gmc, the step-covariance matrix g [equation (|5.17|) ] and gMC, an d the 
normalized correlation matrix between steps s(t) and the deviation from the average 
population e(t) — n(t) — n*, C [equation (|5.22() ] and Cmc- All numbers are given to 
four significant digits. Results for TV = 4 are shown in table [3 



TV 



3 



M 



0.9448 
0.8563 



0.7497 0.9450 
0.8935 0.6212 
0.6474 0.9881 



n*/N 



0.8119 
0.7359 



0.6062 
0.4897 
0.5388 



(n)Mc/N 



0.8112 
0.7350 



0.6055 
0.4892 
0.5378 



G 



3.294 -0.8722 
-0.8722 3.224 



3.667 -0.9323 -0.9558 
-0.9323 3.551 -0.9249 
-0.9558 -0.9249 3.590 



G 



MC 



3.295 -0.8784 
-0.8784 3.227 



3.690 -0.9440 -0.9584 
-0.9440 3.573 -0.9289 
-0.9584 -0.9289 3.592 



4.455 1.368 
1.368 3.882 



3.039 0.7902 0.8768 
0.7902 2.285 0.7152 
0.8768 0.7152 2.592 



gMC 



4.442 1.352 
1.352 3.865 



3.050 0.7865 0.8716 
0.7865 2.287 0.7098 
0.8716 0.7098 2.586 



-2.227 -0.6494 
-0.7189 -1.941 



-1.520 
-0.2649 
-0.5956 



-0.5253 
-1.143 
-0.1908 



-0.2812 
-0.5244 
-1.296 



^MC 



-2.221 -0.6385 
-0.7139 -1.933 
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Table 2. Theoretical results for a TSS community with JV = 4, compared with 
corresponding quantities from a Monte Carlo simulation with No — 10, 000 over 524,290 
generations. The quantities shown are the same as in table ^ The eigenvalues of S 
are relatively close to unity, so evaluation of G from equation (|A.36|) required 180 
iterations. 



M 



/ 0.5507 
0.9543 

0.4724 0.6190 

\ 0.9734 0.6808 



0.5101 0.9847 \ 

0.8437 0.9508 

0.7371 
0.1862 



n*/N 



/ 0.3355 \ 
0.7034 
0.2366 
V 0.3468 / 



(n)Mc/N 



/ 0.3356 \ 
0.7035 
0.2354 
V 0.3472 J 



G 



/ 4.215 -1.043 -2.583 1.132 \ 

-1.043 3.882 -0.1463 -1.006 

-2.583 -0.1463 6.509 -3.791 

V 1.132 -1.006 -3.791 5.653 / 



G 



MC 



/ 4.252 
-1.055 
-2.620 

V 1.151 



-1.055 -2.620 1.152 \ 

3.870 -0.1366 -1.003 

-0.1366 6.612 -3.879 

-1.003 -3.879 5.718 



/ 1.387 0.6468 0.2384 0.3094 \ 

0.6468 3.704 0.4482 0.6518 

0.2384 0.4482 0.8927 0.2414 

\ 0.3094 0.6518 0.2414 1.461 / 



gMC 



/ 1.388 0.6417 

0.6417 3.694 

0.2368 0.4440 

V 0.3090 0.6501 



0.2368 0.3090 \ 

0.4440 0.6501 

0.8870 0.2381 

0.2381 1.462 / 



/ -0.6934 
-0.2517 
-0.1655 

V -0.1253 



-0.3951 
-1.852 
-0.2681 
-0.3237 



-0.07285 
-0.1801 
-0.4464 
-0.2209 



-0.1841 \ 
-0.3281 
-0.02048 
-0.7304 / 



/ -0.6942 
-0.2510 



V 



-0.1642 
-0.1239 



-0.3907 
-1.847 
-0.2661 
-0.3242 



-0.07252 
-0.1779 
-0.4435 
-0.2202 



-0.1851 \ 
-0.3259 
-0.01783 
-0.7310 / 
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and the correlation of these steps with the deviations from the average: e(t). 

In a stationary state, the mean-field prediction for the step size is clearly zero 
(s MF = 0), as must also be the case for the average (s)* . Nevertheless, we expect a 
typical step size to be O (V^o) • For a more detailed picture, we may seek the step-size 
distribution (SSD) in the steady state: V* (s). A precise definition is 

V* (s) = lim V 5 (s, ft -n)V (ft, t + 1; n, t) , (5.2) 



t—>CO ■ 

n'n 



where V (ft, t + 1; ft, t) is the joint probability for finding the system with population 
ft at time t and with ft in the next step. From the master equation f)H.7j) . it is clear that 
this is just R (n'\n) V (ft, t), so that 

UmV(n',t + l;ft,t) = R(ft\ft)V* (n) . (5.3) 

t— >oo 

Thus, once the steady-state distribution V* (n) is known, the SSD can be obtained from: 
■p* (s) = J2R(n + s\n) V* in) . (5.4) 

n 

Applying this formalism to our case of coevolving species without mutations, we 
can exploit all the approximations detailed in the previous section, namely, a Gaussian 
for the stationary state (V* (ft) = Vq (n) of equation (|4.5|l ). continuous variables for 
n, and integrals instead of sums. Given equation ()4.2j) . the success of that scheme is 
implicitly dependent on the fact that R(n'\n) is also well approximated by a Gaussian 
|| . As a result, we need not carry out another lengthy analysis to conclude that V* (s) 
should also be of the form 

1 



VI (s) = (2nN y Af/2 (det 7) 1/2 exp 



2N< 



■sa IJ sj 



(5.5) 



However, it is clear that the matrix 7 is distinct from f , since the former must contain 
some "dynamic" information. (Note that the caret is to remind us that, despite the 
presence of two superscripts, j IJ is of 0(1).) Now, within the context of simple 
Gaussians, 7 can be found by computing its inverse 

g = 7 _1 (5.6) 

which is just the second moment 

s 

Let us first derive an exact formula for this quantity in terms of (•)*. Starting with 

s,n 

= K - n i) i n 'j - nj) R {n'\n) V* {ft) , (5.9) 

ft' ,ri 

|| In the same manner as for the stationary distributions, this property can be derived from the definition 
of R using straightforward, but tedious, manipulations with the binomials. Here, we can treat this as 
an assumption, the justification of which will be the agreement with simulation data. 
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we use equations (|4.2j) and ([3.11)1 and arrive at the exact relation: 

(s lSj ) s = 2 (njnj)* - (fmj)* - (njfj)* . (5.10) 

Next, we exploit V* (n) = Vq (n) for computing (•}* and apply equation (|A.17J) . 
Thus, 

(n inj )* = nifij + N G U (5.11) 

and 

(finj)* = fj (n) nj + N fj K (n) G KJ + ^f™ (n) njG KM ■ (5.12) 

Thanks to equation (jA.18|) . the first and the last terms in equation (|5.12j) combine, so 
that 

(fmj)* = funj + jVo/f (n) G KJ . (5.13) 

To the order kept here, (n) is just ff (n*) = Sf. So, collecting various items and 
using equation ([3.25)1 . we arrive at 

(sisj) s S N (2Gjj - SfG KJ - G IK S«) (5.14) 

= -No {k?G K j + G IK Af^j . (5.15) 

Finally, we relate this result to the Gaussian approximation equation ()5.5j) intended for 
the SSD, which provides 

(sisj) s = N o9lJ , (5.16) 

and obtain a simple equation for g: 

g = -(AG + GA T ) . (5.17) 

As in the previous section, these predictions are well borne out in simulations. For 
the same M = 2 above, we display a two-dimensional histogram of the step 

sizes in figure El Note that this distribution is indeed quite different from that for the 
populations (figure [TJ). For the Af > 2 cases, full histograms are difficult to display and 
we only show the correlation matrices g and their Monte Carlo counterparts gMc m 
tables □ and 

Although the SSD probes the underlying dynamics, it does not contain all 
information of the stochastic (Ornstein-Uhlenbeck) process. In particular, equation 
(I5.17J) shows that the antisymmetric part of AG is "missing." To remedy this 
shortcoming, we turn to another question which naturally comes to mind: How are the 
steps (s) correlated with the deviations from the average (e = n — n)7 Note that, since 
two quantities are involved (s and e), a full distribution associated with this question 
is slightly more involved than V* (n) or V* (s). Nevertheless, within the approximation 
scheme we use, it would be a (generalized) Gaussian in the steady state. For simplicity, 
let us focus on, instead of the full distribution, only the normalized correlation of s with 
e: 

C u = ( Sl sj) /N . (5.18) 
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Figure 2. Contour plot of the theoretical Gaussian and simulated joint probability 
densities for the steps, si and s%, equation (|5.5() . The line types and parameters are 
the same as in figure ^ an d the simulated histogram was obtained from the same 
simulation run as in that figure. The positive slope of the long axis indicates that the 
steps are positively correlated. 



By using this notation, we are displaying again our expectation that {sjSj) is also 
of O (N ). Starting at the same point as equation (|5.8|) . we write 

(s l£j ) = ^2s I £ J R(n + s\n)V* (ft) (5.19) 

= K - ij) (nj - nj) R (n\n) V* (n) . (5.20) 

Since (ft')* = (ft)*, this reduces to 

(siej) = (fmj)* - (n/nj)* , (5.21) 

which is again an exact relationship. Repeating the analysis in the previous paragraph, 
we find, in the Gaussian approximation, a simple result: (siEj) = N AfG KJ , or 

C = AG. (5.22) 

Unlike the step-size covariance, this matrix is not symmetric in general and contains the 
full information of the dynamics (to this order of our approximation). Note also that, 
since A is negative and G is positive, this correlation is typically negative. This simply 
reflects a restoring dynamics: steps and deviations from the mean tend to be opposite. 
In tables ^ and 121 we see that there is excellent agreement between the predicted C's 
and their counterparts from simulations. 

Another commonly studied correlation is the connected two-time function, 
({nj(t + r)nj(t)) — (n/)(nj)) /N , which, in the stationary state, is just 

E u {r) = (si{t + r) £l {t))/N . (5.23) 
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Within the Gaussian approximation, equation (|4.18|) can be used repeatedly to express 
e(t+r) in terms of e(t) and the noise, rjj(t + t') at each time step. Since { r q I {t + t')e{t)) = 
{r]i{t + t'))(e{t)) = for all t', we find 



6. Concluding remarks 

Since we have presented a considerable amount of mathematical details above, it is 
worthwhile to provide a short summary. In the model presented in jl], the entire 
evolution depends only on the parameters, fj,, F , N , and Mjj. Strictly speaking, 
the only stationary state corresponds to total extinction (n = 0). But, if mutation 
is suppressed (/i = 0), then, not only does the configuration space break up into 
many sectors, but there will also be non-trivial long-lived (O (e^ ) generations) states. 
Referred to as "truly stationary states," these communities consist of a fixed number 
(Af > 1) of species. Typically, the population of each species is O (N ). Provided we are 
not near a "critical point," the full distribution of these populations can be studied by a 
systematic approximation scheme, starting with a multivariate Gaussian, parametrized 
by its mean and covariance matrix. For each TSS, our theory predicts these parameters 
and thus, the full distribution. Specifically, from the set of parameters, fi, F , N , 
and Mjj, we can compute ft*, H, and S = 1 — A ( equations 13.23113. 26|) . Then, at the 
lowest order in our approximation scheme, the mean n and the covariance matrix G 
are given explicitly (equations 14. 101 and 14. 16114. 17J) . In addition to this "static" aspect 
of the steady state, we also presented two ways of characterizing the "dynamic" aspect. 
One is the full distribution of sizes of steps (changes in the populations in a single time 
step, denoted by s). Well approximated by a Gaussian also, this distribution has zero 
mean and covariance g , which is given explicitly by equation (|5.17|) . The other is 
the correlation of the steps with the populations just before the step (specified by the 
deviations from the average, e). Denoted by C, this correlation is given explicitly by 
equation (|5.22|) . Finally, we have shown that there is excellent agreement between these 
analytical predictions and Monte Carlo simulations in ten typical TSS communities 
(data for three of which are shown in tables Q and El and figures Q and EJ) • 

With a solid understanding of the steady-state properties of stable communities in 
the absence of mutations, our next goal is a study of quasi-steady states of populations 
with low mutation rates. To carry out a serious analysis with \i > 0, we must define 
the restrictions on R carefully. Otherwise, it would be impossible to find a satisfactory 
solution to the equation 



As in the case for fi = 0, we may seek an approximate solution, in the form of a product 
of a Gaussian distribution for the dominant species and exponential distributions for 
the minority mutants [16J. Provided this program is successful, the next step would be 
to study the probability of this type of QSS having "dangerous" mutants at the higher 



3(r) = S r G . 



(5.24) 




(6.1) 
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order of /i, which hopefully will lead to some understanding of the distribution of QSS 
lifetimes [12] and the presence of 1/ f noise in power spectral densities. Beyond this step, 
perhaps sophisticated renormalization-group techniques can be marshalled to account 
for the self-similar patterns displayed, as well as to identify universality classes for such 
behaviour. Needless to say, even for such a simple model of coevolution, much remains 
to be explored. 
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Appendix A. Mathematical Detail 

Starting from the general expression for a discrete-time Markov process, 

V (ft, t + 1) = R (#1") ? & *) , (A.l) 

n 

we write the mean-field approximation for the evolution equation as 

n MF (t + l) = f(n MF (t)) , (A.2) 

where the functional form of / is given by 

f(n) = ^n'R(n\n) . (A.3) 

n> 

Note that the components of ft play the role of coordinates in equation (|A.1|) and take 
only non-negative integer values in models of population dynamics. However, there is 
no guarantee that / will be an integer in equation (|A.2|) . so that we must allow n MF 
to be continuous variables (functions of t). As we will see below, our analysis is much 
simplified if we also assume n to be continuous. 

Focusing only on simple fixed points (as opposed to fixed cycles involving two or 
more points) of the mean- field theory, we denote a FP by ft*. (To keep the notation 
from being too cumbersome, we suppress the superscript MF and only keep in mind that 
n* represents just M real numbers.) The equation for determining n* is 

n* = f (ft) . (A.4) 

Anticipating the need for the associated stability matrix, we write it as 

S = Vf , (A.5) 
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where V stands for derivatives with respect to ft . Explicitly, the matrix elements of S 
are 

Sj = dfj/dnj^ , (A.6) 

clearly asymmetric in general. As in our earlier study, we restrict ourselves to stable 
fixed points, so that the real parts of the eigenvalues of S (denoted by 1 + A in jl]) 
should lie in the interval (—1, 1). Indeed, we restrict our attention here to isolated FP's 
far from others, so that no (real part of any) eigenvalue is close to unity. 

Turning to the full stochastic problem, we consider a stationary distribution V* (ft), 
which satisfies 

V* (n') = J2 R ( H '\ H ) V * • ( A - 7 ) 

n 

Therefore, it is a right eigenvector of the matrix R(n'\ft) with unit eigenvalue. Its 
existence is guaranteed (though not necessarily its uniqueness, in general) by the 
presence of a left eigenvector of unit eigenvalue (i.e., u(ft) = 1 Vn), thanks to having 
probability conserving rates. If V* is known, then we can compute the stationary 
averages of all quantities via equation (j4.3|) . A different approach is to study the 
equations satisfied by these averages. Multiplying equation (|A.7|) by a quantity and 
summing, we find, in general, no closed equation. Instead, a tangle of infinitely many 
coupled equations appear. Examples are 

(niY = {f j (ft))* and (njnj)* = (H u (ft))* , (A.8) 

where // is an element of /, given by equation ljA.3|) . and 

H T j (ft) = ^rijrijR (n'\ft) . (A.9) 
n> 

In both cases, the right-hand sides contain expectation values of all powers of ft, so that 
infinitely many equations are needed to close the set. 

In practice, for generic situations such as those described in the main text, we may 
approximate V* by a Gaussian: 



1 /2 

V* (ft) = V G (ft) = (2rrN y Af/2 (det f J exp 



2Ni 



(m - m) f IJ (rij - nj) 



(A.10) 



where the unknowns nj and T IJ are to be determined. Naturally, we expect that the 
center of the Gaussian lies close to the mean- field FP, ft*. This scenario is corroborated 
by our case studies above and our analysis below. Indeed, this approximation can be 
used as the starting point of a systematic expansion, which relies on having a large 
parameter (No) controlling the fluctuations to O(y/No). In this context, n and T 
are assumed to be O (N ) and 0(1), respectively. In addition to playing the role of 
ordering a systematic expansion, a large Nq (and so, large ft*) provides two further 
simplifications: (i) n may be promoted to continuous variables so that sums can be 
replaced by integrals, and (ii) integration can be extended to [— oo, oo] while only 
incurring errors of O ^e~ v/ ^° 
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Once fij and T IJ are determined, then we can set up expansions for the averages 
of any quantity Q: 

(Q(n))* = Y^Q{n)V* [ft) , (A.ll) 

n 

starting with 

Q{n)V G [ft) dn . (A. 12) 



Given that nj does not deviate far from ni, such integrals can be handled by first 
expanding Q around nf. 

Q{n) =Q + Q I (n I - n,) + X -Q 1J (m - n T ) (nj - nj) + . . . , (A. 13) 

where the summation convention is used and 

Q = Q{n)\ a=n , (A. 14) 

dQ {ft) 



Q 1 
Q 



dni 
u_ d 2 Q(n) 



(A.15) 
(A.16) 



dnjdnj 

The integral can now be performed, so that 

/AT 
Q(n)V G (ft) dn = Q + -fQ IJ G u + ... , (A.17) 

where G = f 1-1 . There is no need to be alarmed at the factor iVg in the second term. 
Since Q IJ involves two derivatives with respect to n, it is typically of the order of 1/Nq 
compared to Q. Thus, the second term in this expansion is O (1/JVo) compared to the 
first. 

With this machinery, we seek equations to fix fij and G. From equation (|4.4|) . we 
have, with Q = fi, 

nj = (fj (ft))* = fj (n) + ^fj JK G JK + ... . (A.18) 

Again, we expect the second term to be O (1/N ) compared to the first, so that we may 
write an expansion for fif. 

nj = nf +n i j ) + ... . (A. 19) 

Inserting this back into equation (|A.18|) and noting equation (|A.4|) . we find the expected 

nf = n} . (A. 20) 

Proceeding to the next order, we find 

n? = H [nf) n? + ^fl^JK , (A.21) 

where all functions of n in the last term can be evaluated at nf\ To save notation, we 
will just write Gjk in lieu of the more explicit form: Gj K . We recognize that // 
is just the stability matrix Sj and write 

-(1) = ^l V M fM JK djK (A22) 
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where 

V = (l-S) _1 (A.23) 

(the inverse of —A in jlj). Let us emphasize that V T M and Gjk are supposedly O (1), 
while fu 3K is O (1/Nq) so that the right-hand side of equation (jA.22|) is a quantity of 
0(1). This expression also provides a precise meaning to the phrase "no eigenvalue (of 
S) is close to unity," which appeared as a restriction on the FP's we study. 
Extending this technique to njnj , we consider 

NoGu + (nj)* (nj)* = (mnj)* = (H u (ft))* , (A.24) 

from equations (|4.4|) and IjA.Sft . Again, using the Gaussian approximation for both sides, 
we arrive at 

N Gjj + njnj = H n (n) + ^Hjj KM G km + ... . (A.25) 

At the lowest order, 0(Nq), of this equation, internal consistency will ensure that 
Hjj (n^) = fiyfi^ and should provide a check for tedious, error-prone computations. 
Furthermore, we can consider the difference : 

H u (n) - fi (n) fj (n) (A.26) 

and find that it is often one order lower (as explicitly shown in the main text for our 
model). In other words, we have 

H u (ft*) - // (ft) fj (rf ) ~ O (N ) . (A.27) 

To emphasize this property, we explicitly extract a factor N and define 

Hu (ft) = [H u (ft) - fr (ft) fj (ft)} /N , (A.28) 

so that 

H u (n) = fi (n) fj (n) + N H U (n) . (A.29) 

Since we keep only terms of orders Nq and A'o, we can evaluate Hu (n) at n* and define, 
for simplicity, 

H u = Hu(n*) . (A.30) 

As a reminder, this quantity is of O (1). Next, to lowest order needed (i.e., O (1)), we 
find 

tj KM qK qM i qM qK i r r KM , r KM r i a q i \ 

Hu = bj bj + b f bj + fifj + fi fj , (A.31J 
all evaluated at the FP. Inserting these expressions into equation (|A.25|) . we obtain 
Nod j + run; S fj (n) fj (n) + N Hu 

, Nq r K M M K RM KM -, 

+ [ b I b J + b I b J + III J + // JJ\ ^KM + ■ • • • 

(A.32) 



Exploiting equation (jA.18|) . we arrive at an equation for G : 

Gu = H u + \ [SfSf + SfSf] G KM ■ (A.33) 
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But, G is symmetric, so that the final equation, in matrix form, reads 

G - SGS T = H . (A.34) 

Here, S T denotes the transpose of S, and, from equation (|A.9J) . H is necessarily 
symmetric. Again, let us emphasize that all quantities in this equation are 0(1). 
Unfortunately, there is no simple way to express G in terms of the known matrices: 
S and H. A formal series can be written as 

G = H + SHS T + SSHS T S T + ... , (A.35) 

equivalent to a recursion relation, 

G n = H + SG n _!S T (A.36) 

with G = H. It was this recursion method that was used to obtain numerical results 
for comparison with the Monte Carlo results in the figures and tables. Convergence to 
four significant digits was obtained with n = 20 for M = 2, n = 30 for M = 3, and 
n = 140 for M = 4. 

An explicit form for G can be found when we examine this equation in the frame 
where S is diagonal. To be explicit, take matrix elements of equation (|A.34|) with u^, 
the left eigenvectors of S, which satisfy 

u ! a Sf = Xa^a (no sum on a) . (A. 37) 



The result is 



where 



G ab = Habl (1 - A a A b ) (no sum) , (A.38) 



G ab = Guu^ui and H ab = Huu^ui . (A.39) 

To find the original matrix elements, we apply the dual set vf (i.e., the right eigenvectors 
of S, normalized by v^u{ = 5%) and obtain 

G u = v*v b jG ab = ^Habl (1 - A a A b ) . (A.40) 

a, b 

Finally, this expression again shows the importance of insisting that "no eigenvalue 
(of S) is close to unity." In case one or more A approaches unity, the system would be 
labelled "critical" or "multi-critical," in the language of phase transitions. Away from 
such points, the Gaussian approximation, along with its associated Ornstein-Uhlenbeck 
process, is adequate. For further details, see, for example, the books by van Kampen ^H] 
and Risken [T7] . 
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